【Python金融量化】VaR系列（四）：蒙特卡洛方法估计VaR

作者：量化小白H     Python爱好者社区专栏作者
个人公众号：量化小白上分记

前文传送门：

之前几篇总结的方法都是对于向前一日VaR的建模，都以是以VaR=波动率乘以分布函数逆函数为基础。但如果要计算向前k日的VaR，如果还使用上述公式，波动率和分布函数应该换成k日滚动窗口的，好像还没见过这样的Garch模型。

所以本文总结一种可以对向前k日VaR进行估计的方法，蒙特卡洛方法。需要说明的是，通过蒙特卡洛方法估计VaR可以从很多角度入手，本文只是其中一个角度，并非唯一建模方法。

01

模型推导

Monte-Garch

前面提到，之前的基础公式不太适用K日VaR的估计，需要寻找其他方法。蒙特卡洛方法从VaR定义出发，要计算资产在给定概率下，未来K日的最大损失，我们可以模拟出资产未来k日的日收益率，计算得到资产的K日收益模拟序列，多次模拟后，类似之前HS和WHS方法，取整个模拟序列的p分位数即可，具体操作如下

我们依然假设资产的日收益率序列服从正态分布，从而标准化收益率服从标准正态分布。

同时收益率序列的波动率可以用Garch族模型进行建模，假设波动率服从Garch(1,1)模型

这样，我们可以用标准正态随机数对向前一天的标准化收益率z进行模拟，记模拟次数为MC。

从而可以得到模拟的收益率为

同理对向前两天的标准化收益率z进行模拟MC次

波动率通过Garch模型进行更新

之后可以得到向前两天的日收益率

如果要估计向前K天的VaR，我们将上述过程重复K次，过程可以表述如下

这样向前K日的总收益率可以把k日收益率加总得到（如果是算单利的话，复利用乘法）

这样我们就得到了对于未来K日总收益的MC个模拟 值，我们用这MC个模拟值估计未来K日的VaR，类似HS方法，直接取p分位数即可

Monte-FHS

上面的方法是将蒙特卡洛模型与Garch组模型相结合，此外，也可以将蒙特卡洛方法与历史模拟方法相结合，即不对未来收益率进行建模，而是直接用过去一段时间的收益率作为未来每日收益率的估计量，对上述方法做一点小的修改，具体如下

对于过去一段时间的日收益率，我们可以通过Garch模型对波动率进行估计，得到过去一段时间的标准化收益率

随后，我们每次都从这个收益率序列中进行有放回抽样，生成对未来k日的收益率模拟序列，再重复波动率的估计和收益率的加总过程

最后求出VaR

这里做一点分析，这两种方法可以认为属于两类方法：一类可以认为是完全没有模型的方法，包括WHS,HS，FHS，对于收益率的分布不做任何估计，遵循历史会重演的原则，用历史收益率作为未来收益率的估计量，一类是有严密假设前提和模型推导的参数方法，包括Garch，EVT等等。

两种方法各有优劣，从之前文章的结果来看，没有模型的方法基本上是不能看的，不过从本文之后的两种方法的对比来看还可以。而有假设有模型的方法结果可能会比较好，但前提是能对收益率序列作出合理假设，从正态分布，t分布，渐近t分布等等，理论研究中尝试了各种分布，但显然没有一定最有效的分布，一切都在于尝试。

Monte-DCC-Garch

刚才两种方法都是对单个资产的VaR估计，也可以把蒙特卡洛方法与前一篇文章中的DCC方法相结合，估计组合的向前k日VaR。用Monte-Garch或者Monte-FHS都可以，过程差不多，这里以Monte-Garch为例。

组合VaR估计与单资产VaR估计的不同在于组合不仅需要估计资产的波动率，还需要估计资产之间的相关性，换句话说，需要估计资产的协方差阵。和之前类似，我们可以模拟组合的日标准化收益率序列，用Garch模型更新波动率，用DCC模型跟新相关系数，然后计算组合的日收益率，K日总收益，最后计算组合的向前K日VaR。

但如果我们直接用之前的方法分别模拟每个但资产的收益率，显然各个资产之间是不相关的，因此我们必须得到符合给定相关系数矩阵的模拟收益率序列，以2个资产为例进行说明

根据之前的方法我们可以得到如下的两个不相关的标准化收益率序列

我们的目标是得到符合给定相关系数rho的标准化收益率序列

显然通过如下转换即可完成

两资产情况下，相关系数矩阵的平方根很容易算出来

多资产情况下，相关系数矩阵的平方根不是很好算，但也可以算出来，用到cholesky分解，不过如果是写代码的话，也就是一行代码调个函数的事情。

综上，计算向前K日VaR整个过程可以表述如下：

1. 模拟生成不相关的资产标准化收益率序列，均服从标准正态分布

2. 将不相关的资产标准化收益率序列转换为相关系数矩阵为给定值的标准化收益率序列

3. 将标准化收益率转化为收益率，按照各资产的权重进行累加，得到组合的收益率

4. 通过Garch模型更新波动率，通过DCC模型跟新相关系数矩阵

5. 重复1-4过程K次，将K次的收益率序列相加总收益率

6. 重复1-5过程MC次，得到MC个总收益率，取p分位数即可。

02

实证分析

软件：python3.6

数据：S&P500、US 10yr T-Note Fixed Term（同上一篇）

区间：2001-2010

蒙特卡洛模拟次数：10000次

数据和代码在后台回复“VaR4”获取

单资产的VaR估计均在数据序列最后一天向前估计，估计K=500天

波动率初始值使用当天的波动率，波动率建模用之前文章中的Ngarch（1，1）模型，两资产组合的VaR估计类似

Monte-Garch（代码见VaR_Monte_NGarch）

上一篇中提到，Garch模型具有均值回归的特性，最终会回到长期平均水平，如果我们将波动率初始值分别改为高初始波动率：3倍长期平均波动率和低初始波动率：0.5倍长期平均波动，结果如下

高初始值的结果先升高再降低，低初始值的结果持续升高，这个结果看起来不是特别明显，文献中的文献中的结果要比我的好很多

Monte-FHS（代码VaR_Monte_FHS）

对比Monte-Garch和Monte-FHS的结果，整体趋势差不多，不过Monte-FHS得到的值要相对小一些。不同初始波动率的结果对比如下

此外，代码中同时还对比了用RM，Monte-Garch、Monte-FHS计算出的向前一日VaR，其中

RM：0.045864978982555343；

Monte-Garch：0.052665510343511503；

Monte-FHS：0.053369497694956955

组合VaR估计（代码见getVaR_Monte）

这里与上一篇中DCC模型估计的向前一日VaR结果进行对比

黑色为Monte-Garch-DCC的结果，红色为DCC的结果，几乎是完全重合的，黑色被覆盖。但效率上讲，蒙特卡洛的运行速度较慢，远不如DCC。

03

代码

这里只给出主要函数代码，全部代码后台获取，另外代码实现的过程实际是参考文献中Chapter8中的习题，注释写的不多，看不明白的可以翻一翻文献，如果想把代码用到别处，需要做一些修改，比如组合收益率中，因为这里考虑的是等权重组合，所以并没有设置权重参数，直接取的平均，如果要考虑非等权，函数输入参数中需加上权重参数。

• Monte-Garch

 1def VaR_Monte_NGarch(num,p,sigma_first,ndays,omega,alpha,beta,theta):
 2    np.random.seed(4)
 3    MC = num
 4    sigma2_Garch = sigma_first
 5    data_MC = pd.DataFrame(index = range(MC))
 6    data_R = pd.DataFrame(index = range(MC))
 7    data_MC['z_day1'] = np.random.normal(size = (10000,1))
 8    data_R['R_day1'] = sigma2_Garch**0.5 * data_MC['z_day1']
 9    data_MC['sigma2_day1'] = omega + alpha*(data_R['R_day1'] - theta*sigma2_Garch**0.5)**2  + beta*sigma2_Garch
10
11    # 低2-10天公式一样，循环
12    for i in range(2,ndays +1):
13        exec("data_MC['z_day" + str(i) + "'] = np.random.normal(size = (10000,1))")
14        exec("data_R['R_day" + str(i) + "'] = data_MC['sigma2_day" + str(i-1) + "']**0.5 * data_MC['z_day" + str(i) + "']")
15        exec("data_MC['sigma2_day" + str(i) + "'] = omega + alpha*(data_R['R_day" + str(i) + "'] - theta*data_MC['sigma2_day" + \             str(i -1) + "']**0.5)**2 + beta*data_MC['sigma2_day" + str(i-1) + "']")
16    VaR = pd.DataFrame(index = range(ndays))
17    VaR['ndays'] = np.arange(1,ndays+1)
18    VaR['VaR'] = 0
19    for i in range(ndays):
20        R_ndays = data_R.iloc[:,:i+1].sum(axis = 1)
21        VaR.loc[i,'VaR'] = - np.percentile(R_ndays,p) 
22    return VaR

• Monte-FHS

 1def VaR_Monte_FHS(num,p,sigma_first,ndays,omega,alpha,beta,theta):
 2    np.random.seed(4)
 3    MC = num
 4    sigma2_FHS = sigma_first
 5    data_FHS = pd.DataFrame(index = range(MC))
 6    data_FHS_R = pd.DataFrame(index = range(MC))
 7    data_FHS['z_day1'] = data.loc[np.random.randint(0,data.shape[0]-1,MC),'z_Garch'].values
 8    data_FHS_R['R_day1'] = sigma2_Garch**0.5 * data_FHS['z_day1']
 9    data_FHS['sigma2_day1'] = omega + alpha*(data_FHS_R['R_day1'] - theta*sigma2_FHS**0.5)**2  + beta*sigma2_FHS
10
11    # 低2-10天公式一样，循环
12    for i in range(2,ndays +1):
13        exec("data_FHS['z_day" + str(i) + "'] = data.loc[np.random.randint(0,data.shape[0]-1,MC),'z_Garch'].values")
14        exec("data_FHS_R['R_day" + str(i) + "'] = data_FHS['sigma2_day" + str(i-1) + "']**0.5 * data_FHS['z_day" + str(i) + "']")
15        exec("data_FHS['sigma2_day" + str(i) + "'] = omega + alpha*(data_FHS_R['R_day" + str(i) + "'] - theta*data_FHS['sigma2_day" + \             str(i -1) + "']**0.5)**2 + beta*data_FHS['sigma2_day" + str(i-1) + "']")
16
17    VaR = pd.DataFrame(index = range(ndays))
18    VaR['ndays'] = np.arange(1,ndays+1)
19    VaR['VaR'] = 0
20    for i in range(ndays):
21        R_ndays = data_FHS_R.iloc[:,:i+1].sum(axis = 1)
22        VaR.loc[i,'VaR'] = - np.percentile(R_ndays,p) 
23    return VaR

• Monte-DCC-Garch

 1def getVaR_Monte(n,alpha_DCC,beta_DCC,sigma_first,rho_first,alpha_garch,beta_garch,omega_garch,theta_garch,p):
 2    MC = 10000
 3    sigma2_Garch = np.dot(np.diag(sigma_first),np.ones([2,MC]))
 4
 5    np.random.seed(4)   
 6
 7    VaR = pd.DataFrame(index = range(n))
 8    VaR['ndays'] = np.arange(1,n+1)
 9    VaR['VaR'] = 0   
10    VaR['rho'] = 0
11    VaR.loc[0,'rho'] = rho_first
12
13
14    # 记录组合的模拟收益
15    data_R = pd.DataFrame(index = range(MC))
16
17    # 模拟组合中各资产的标准化收益率
18    z = np.random.normal(0,1,size = (2,MC))
19    gamma_half = np.array([[1,0],[VaR.loc[0,'rho'],(1-VaR.loc[0,'rho']**2)**0.5]])
20    # 转换为指定相关系数的序列
21    z = np.dot(gamma_half,z)
22
23    # 各资产的真实收益率 = 标准化收益率*波动率
24    r = np.multiply(sigma2_Garch**0.5,z)
25
26    # 两资产等权重，总收益率
27    data_R['R_day1'] =  np.mean(r,axis = 0)
28
29    if n==1:
30
31        R_ndays = data_R.iloc[:,:2].sum(axis = 1)
32        VaR.loc[0,'VaR'] = - np.percentile(R_ndays,p) 
33        return VaR
34
35    # 若计算向前多天的VaR,通过DCC,GARCH模型更新参数
36
37    alpha_garch = np.dot(np.diag(alpha_garch),np.ones([2,MC]))
38    beta_garch = np.dot(np.diag(beta_garch),np.ones([2,MC]))
39    omega_garch = np.dot(np.diag(omega_garch),np.ones([2,MC]))
40    theta_garch = np.dot(np.diag(theta_garch),np.ones([2,MC]))    
41    # 波动率更新 - Garch模型
42    sigma2_Garch = omega_garch + np.multiply(alpha_garch,(r - np.multiply(theta_garch,sigma2_Garch**0.5)))**2 + np.multiply(beta_garch,sigma2_Garch)
43
44    # 相关系数更新 - DCC模型
45    q11 = pd.DataFrame(index = range(MC))
46    q12 = pd.DataFrame(index = range(MC))
47    q22 = pd.DataFrame(index = range(MC))
48    q11['day1'] = 1 + alpha_DCC*((z[0,]**2)-1)  
49    q12['day1'] = VaR.loc[0,'rho'] + alpha_DCC*(z[0,]*z[1,]- VaR.loc[0,'rho']) + beta_DCC*(VaR.loc[0,'rho'] - VaR.loc[0,'rho'])
50    q22['day1']= 1 + alpha_DCC*((z[1,]**2)-1)
51
52    VaR.loc[1,'rho'] = np.mean(q12['day1']/((q11['day1']*q22['day1'])**0.5))
53    # gamma更新
54    gamma_half = np.array([[1,0],[VaR.loc[1,'rho'],(1-VaR.loc[1,'rho']**2)**0.5]])
55
56
57    for i in range(2,n):
58        # 生成不相关的随机序列，并转换为指定相关系数的随机数
59        z = np.random.normal(0,1,size = (2,MC))
60        z = np.dot(gamma_half,z)
61        r = np.multiply(sigma2_Garch**0.5,z)
62        data_R['R_day' + str(i)] =  np.mean(r,axis = 0)        
63
64        sigma2_Garch = omega_garch + np.multiply(alpha_garch,(r - np.multiply(theta_garch,sigma2_Garch**0.5)))**2 + np.multiply(beta_garch,sigma2_Garch)
65
66
67        q11['day'+str(i)] = 1 + alpha_DCC*((z[0,]**2)-1) + beta_DCC*(q11['day'+str(i-1)]-1)
68        q12['day'+str(i)] = VaR.loc[0,'rho'] + alpha_DCC*(z[0,]*z[1,]- VaR.loc[0,'rho'])+beta_DCC*(q12['day' + str(i-1)] - VaR.loc[0,'rho'])
69        q22['day'+str(i)] = 1 + alpha_DCC*((z[1,]**2)-1) + beta_DCC*(q22['day'+str(i-1)]-1)
70
71        VaR.loc[i,'rho'] = np.mean(q12['day'+str(i)]/((q11['day'+str(i)]*q22['day'+str(i)])**0.5))
72        gamma_half = np.array([[1,0],[VaR.loc[i,'rho'],(1-VaR.loc[i,'rho']**2)**0.5]])
73
74    for i in range(n):
75        R_ndays = data_R.iloc[:,:i+1].sum(axis = 1)
76        VaR.loc[i,'VaR'] = - np.percentile(R_ndays,p)                
77
78    return(VaR)

04

总结

写的有点多，总共三个模型，Monte-Garch，Monte-FHS，Monte-DCC-Garch，前两个模型的实证结果与文献结果一致，大可放心，组合的实证文献中没有给出实证结果，是自己算出来的，不保证正确性，不过仅看一日的是和DCC一致的，所以大概率没什么问题，如果觉得有问题，请和我联系，谢谢。

此外还想说明的一点是，计算未来K日的VaR，K如果取的过大，结果没有多大实际意义，越往后越趋近于平稳，因为在模型迭代过程中，每一天的收益率都是模拟的误差不会变化太多，但是波动率和相关系数的误差会随着预测天数增加逐步增大，所以还是关注比较近的几天比较有意义，比如，可以看一下K=500时，DCC模型得到的相关系数变化情况，不同曲线是分别取不同初始值时的结果，越往后预测相关性越小。

最后，本文所有的模型依然是在正态性假设下，没有考虑非正态情况，非正态的情况下需要用到copula模型，在下一篇文章中给出。

参考文献

Christoffersen, Peter F. Elements of financial risk management. Academic Press, 2011.

Python爱好者社区历史文章大合集：

福利：文末扫码立刻关注公众号，“Python爱好者社区”，开始学习Python课程：

关注后在公众号内回复“课程”即可获取：

小编的Python入门免费视频课程！！！

【最新免费微课】小编的Python快速上手matplotlib可视化库！！！

崔老师爬虫实战案例免费学习视频。

陈老师数据分析报告制作免费学习视频。

玩转大数据分析！Spark2.X+Python 精华实战课程免费学习视频。